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Abstract 

We present a practical method to generate classical trajectories with fixed ini- 
tial and final boundary conditions. Our method is based on the minimization 
of a suitably defined discretized action. The method finds its most natural 
application in the study of rare events. Its capabilities are illustrated by non- 
trivial examples. The algorithm lends itself to straightforward parallelization, 
and when combined with ah initio molecular dynamics (MD) it promises to 
offer a powerful tool for the study of chemical reactions. 

The dynamics of complex systems is often characterized by the occurrence of rare but im- 
portant events. Examples of these crucial events are chemical reactions, diffusion processes 
and molecular conformational changes. The reason for the slowness of these processes is 
often related to the existence of large energetic barriers that the system has to clear in order 
to pass from the initial to the final state. Since this transition rate depends exponentially 
on the barrier height, the reaction characteristic time can exceed present-day computational 
capabilities by several orders of magnitude for barriers larger than ~ ksT. Most meth- 
ods are based on the search for the transition state starting from the initial configuration. 
The dynamical properties are then reconstructed by making use of transition state theory. 
In this class of methods one can include the traditional quantum chemical approach [|TJ, 
conformational flooding |2j], hyper-dynamics and the parallel replica method [[|. These 
methods fail when transition state theory is not applicable. Furthermore it often proves 
rather difficult to locate the transition state. This is especially true in complex systems in 
which the potential energy surface (PES) is very rugged and exhibits very many stationary 
points. In such a case the very concept of transition state as a saddle point in the PES is 
called into question. 

Elber and Karplus || have suggested an innovative approach. Rather than starting from 
the initial state and looking for a transition state, they consider a path that connects the 
initial and the final state. This approach is in spirit similar to the more modern nudged 
elastic band approach of Jonsson ||. The method is appealing since it focuses on the global 
character of the path rather than on local properties of the PES. However, the paths do 
not reflect the real dynamics as the system crosses the barrier. Following a seminal work 
of Pratt 0, the importance of determining the real dynamical path has been emphasized 
very recently by Chandler and collaborators || who have proposed a method for statistically 
sampling dynamical paths. This new method is a very significant step towards the study 
of rare events in complex systems, as it requires neither the assumption of transition state 
theory, nor the existence of a single well defined transition state. Its only prerequisite is the 



1 



knowledge of at least one reactive dynamical trajectory. However, determining the initial 
trajectory is often difficult. A general strategy suggested by Chandler and collaborators has 
been to anneal very unlikely trajectories. This procedure is computationally demanding and 
does not guarantee that the system ends up in the desired final state. 

Our work addresses this very general problem: given an initial and a final configuration, 
what are the dynamical paths that connect them? In mathematical terms this is a two- 
point boundary value problem. The standard boundary conditions for Newton's equations 
fix instead the initial values of positions and velocities. 

In principle the way to determine the dynamical path from configuration A to configura- 
tion B has been known since 1744 when Maupertuis proposed the principle of least action, 
published a few years later ||. A more precise mathematical formulation is due to Hamilton. 
In modern terms Hamilton's principle can be written as follows. Given a classical dynamical 
system described by the set of coordinates q, its trajectory q(£) with boundary conditions 
q(0)= q^ and q(r) = q# is determined by locating the stationary point of the action: 

S= f L(q(t),q(i))dt, (1) 
Jo 

where L is the Lagrangean L = T — V and T and V are the kinetic and potential energy 
respectively. By varying S under the condition that the initial and final points of the 



trajectory are fixed [ 10 1 Newton's equations of motion follow. We want to turn this well- 



known principle into a useful computational tool. To this effect, following Gillilan and 



Wilson [[I jjj we discretize the trajectory into P equispaced intervals and estimate the action 

as 

5=EA(i(3L^i) 3 -V'(« J ,-)), (2) 

where qo = , qp = qp , A = t/P is the time interval in which we discretize r and we are 
considering unitary masses. Eq. (0) is obtained by using the most primitive quadrature for 
the action integral and by estimating the velocities from the difference between successive 
points in the trajectory [12] . 



The stationary point of this action is expressed by a set of P — 1 linear equations: 



2dV(q j ) 



2q, - q,_! - q,- +1 - A = 0, (3) 

that is, the stationary point of the discretized action is a Verlet trajectory, that is a trajectory 
that is identical to what would have been generated by the well known Verlet algorithm. A 
straightforward search for the stationary point of S (Eq. (0)) is however very difficult, since 
the action is bounded neither from above nor from below and we are looking for a stationary 
point, not necessarily a maximum or a minimum. Furthermore, the nature of the stationary 
points depends on A. This makes the search for the solution of Eq. (|3]) difficult, since it 
becomes a root-finding problem, which in order to be solved requires a good guess of the 



solution and a Newton- Raphson procedure jl3| . The latter step involves the calculation of 
the second-order derivatives of the potential. 

Alternatively, in order to find the solution of Eq. (||) one can look for the minimum of 
the function: 
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= £ (2q, - q,_, - q J+1 - A 2 ^M) 2 , (4) 

which is very close to the discretized version of the Onsager-Machlup (OM) action, as intro- 
duced by Olender and Elber [TJ|. Locating the minimum of this action again involves the 



evaluation of second-order derivatives. This, for an ab-initio molecular dynamics (MD) or for 
the study of large systems within ordinary MD, is a prohibitively expensive task. Further- 
more, in the practice the trajectories thus generated may have a poor energy conservation, 
especially in the region of the transition state, as we shall show in the following. 

We present here a method which generates accurate trajectories and requires only the 
calculation of the forces. To this effect we shall use the fact that the Newton trajectory has 
to conserve total energy. This is also approximately true for the Verlet trajectories. If A 
is sufficiently small the Verlet trajectories do not strictly conserve energy, but instead the 
energy oscillates around a roughly constant value. We shall take advantage of this property 
of the physically relevant trajectories and supplement the action with a penalty function 
which favors the energy conserving trajectories: 

Q{^E) = S + ^{E J -Ef. (5) 

Here [i determines the strength with which we enforce energy conservation, Ej = (q^ — 
qj + i) 2 / (2A 2 )+V(qj) is the instantaneous energy and E its target value. Our definition of the 
kinetic energy is different from that which is normally used in the Verlet algorithm, in which 
the velocities are estimated by the more accurate central difference method. Nevertheless 
also Ej is approximately conserved and for simplicity's sake we shall use this definition here. 
E corresponds to the yet undetermined physical energy. 

For large \i the second term in eq. (f|) dominates, and the functional 0(q 7 -,£ l ) has 
a minimum. For small /i, on the other hand, S will determine the character of the 
extremum. The two regimes are separated by a critical value fi*. In all the systems studied, 
we have found that there exists a rather large interval of values of // > /i* such that G has 
a minimum and this minimum is globally very close to the Verlet trajectories. 

Whenever the target energy E is not known, we can minimize relative to the q,, and 
to E. This procedure allows us to focus on the physically relevant values of E. 

In order to minimize more efficiently, we follow and make the transformation: 

j A p j A 
qj ■< = QaH (Qb - qu) + Y] a„ sin (7m ), (6) 

T n=l T 

thus automatically satisfying the boundary condition. The advantage of using the a n rather 
than the q, is that the a n have a global character. In practice, we first optimize with 
respect to a relatively small number of a„, thus capturing the global feature of the trajectory, 
and then we add the higher frequencies. Each time we use a standard conjugate gradient 
algorithm to minimize 0. This requires only the evaluation of the forces. The scaling 
is therefore linear in the number of degrees of freedom rather than quadratic as in other 
approaches where second derivatives are needed. 

We illustrate the performance of the functional in a series of examples. The most 



elementary one is a one- dimensional double well potential already studied in ref. [[LJ]. In 
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this simple case it is easy to find a Verlet trajectory oscillating between the two minima. 
If we start from this trajectory and add a random component, we find that for /x = a 
minimum or a maximum cannot be found and the Verlet trajectory corresponds to a saddle 
point. This leads to unstable behavior. We therefore studied the behavior of with respect 
to /x. We found that for /x < /x* (which depends on A) we could not locate a stationary 
point. In particular, from the direct calculation of the Hessian matrix at the stationary 
point, we could extrapolate the change of sign of its lowest eigenvalue exactly at /x*. The 
behavior of Euclidean distance in the close neighborhood of /x* foreshadows the closeness to 
the onset of the instability. 

For larger values of /x, however, we find solutions very similar to the Verlet trajectory. 
In order to measure the accuracy of the © solution, we evaluated the norm of the distance 
(in the (P — 1)— dimensional space of the trajectories) between the Verlet trajectory and the 
minimum, as well as the function O in eq. @. The results shown in Fig. 1 give a comparable 
but not identical estimate for the optimal value of /x. It must be stressed, however, that 
to all interests and purposes, the variational solutions and the original Verlet trajectory are 
globally very close, which makes the precise choice of /x not critical. This is very important, 
since in most cases the exact trajectory is not known; in such cases only the criterion of 
minimal O can be used to optimize /x. 

The second example is a minimization of a trajectory in a two-dimensional configurational 
space, namely in the Mueller potential [15], which was invented as a non-trivial test for 
reaction path algorithms: it has two main minima, and a third minimum that somewhat 
hinders the trajectory from seeing the saddle point. The aim of the optimization is to find 
a physical path between the two deepest minima, possibly locating the saddle point. 

We directly compared our results with those of ref. [14]. Both sets of trajectories lead to 
an apparently satisfactory result. They pass exactly through the saddle point, and the overall 
behavior of the trajectories is clearly physical. Our trajectory has of course a larger value 
of O. However, the O-trajectory shows poorer energy conservation close to the transition 
state, as we see in Fig. 2. This is due to the fact that the optimization of O is global and 
local errors leading to energy non-conservation have little weight. Of course one could also 
supplement O with a penalty function to improve the energy conservation; this would still 
require the calculation of second-order derivatives. 

As a last and far less trivial example we look at a process in which the central atom of 
a seven-atom two-dimensional Lennard- Jones cluster migrates to the surface. This process 



has been studied in detail by Dellago et al. |TB[. Our calculations reproduce their highly 



non-trivial results: in Fig. 3 we show how two of the paths we found correspond to the most 
probable ones found through the elaborate procedure of path sampling in [JITJ. As can be 
seen in Figs. 3 and 4, the dynamical path has to pass several transition states. This is a 
severe challenge for other methods. 

It has been suggested by Zaloj and Elber JET] that the O-trajectories can be used to 
increase the value A beyond what is normal in ordinary MD. The basic assumption is that 
by so doing, the effect of the higher frequencies is integrated out and leads to a Gaussian 
distribution of the errors: 



6j = 2q, - qj _! - q j+1 - A . (7) 

With this in mind, we plot in Fig. 4 the distribution of the errors 6j, defined by (0). As is 



4 



exemplified by the figure, we verified that in the G trajectories, the e» appear to be normally 



distributed. Thus one can implement an accelerated MD scheme in the spirit of ref. |17] . 

In conclusion, we have presented a novel method for generating dynamical trajectories 
with pre-assigned initial and final boundary conditions. The method can be easily imple- 
mented within the Car-Parrinello MD approach, offering a powerful tool for the study of 
chemical reactions, and can be combined with path sampling. Its advantages over previous 
schemes are the fact that it requires only the calculation of the forces, its numerical stability 
and the quality of the trajectories: the latter is a direct product of the energy conservation 
constraint. 

As a side effect, the method can lead to easy localization of the saddle points without 
losing the physical soundness of the solutions. Our approach lends itself to very efficient 
parallelization. Multiple time-scale approaches are natural and feasible. We are sure that 
many other applications of our method will be found and we believe that it is a noteworthy 
advance in the field of molecular dynamics. 

It is a pleasure to thank W. Silvestri, A. Gambirasio, M. Bernasconi, and F. Filippone 
for fruitful discussions and suggestions. The research of D. P. at MPI is sponsored by the 
Alexander von Humboldt Foundation. 
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FIGURES 



FIG. 1. Euclidean distance between the G trajectory and the practically exact Verlet trajectory 
(red, left scale) and OM action value (blue, right scale) as a function of /i in a one-dimensional 
double- well potential of the type U(x) = l/2(x 2 +^4*exp(— ax 2 ) — B), with B = l/a[\og(aA) + 1], 
A = 80 and a = 0.04375 14]. With these parameter values, the two minima are located at ±5.33 
and the barrier height is 14.2. The integration time step is A = 0.31 and the total number of P 
points is 64. The starting guess for the minimization of is a randomized Verlet trajectory. For 
H = /U = 0.39 the Hessian of ceases to be positive definite. In the inset, the Verlet trajectory 
(red) is compared with our solution (blue) at jjL = 100. 



FIG. 2. Upper panel: trajectories (blue crosses) and O-trajectories (red dots) on the Mueller 
potential energy surface [14]. The potential energies of the two lowest minima are -146 and -106, 
and the transition state is located at -41 (same units as in fl4| ). Lower panel: energy as a function 



of time for our trajectory (blue line) and O-trajectory (red line). As in reference [14] we have taken 
a time step of 0.01; the total number of points was 300. The regions of accumulation of blue points 
correspond to the two lowest minima. 



FIG. 3. Two of the paths followed by a two-dimensional cluster of 7 atoms for the migration 
of a central atom to the surface. In Lennard-Jones units the total time is 3 and the time step was 
taken to be A = 0.06. The initial trajectory was a linear interpolation between the initial and final 
points. We show in the picture the sequence of metastable states visited by the system. 

FIG. 4. Potential energy profile (top) and conserved energy (middle) for the path (b) in Fig. 
3. The two intermediate metastable states depicted in Fig. 3 are indicated by arrows. The lower 
panel shows the distribution of the errors €j along the same trajectory. The green line shows a 
gaussian fit. 
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